Scaling Behavior of Threshold Epidemics 
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We study the classic Susceptible-Infected-Recovered (SIR) model for the spread of an infectious 
disease. In this stochastic process, there are two competing mechanism: infection and recovery. Sus- 
ceptible individuals may contract the disease from infected individuals, while infected ones recover 
from the disease at a constant rate and are never infected again. Our focus is the behavior at the 
epidemic threshold where the rates of the infection and recovery processes balance. In the infinite 
population limit, we establish analytically scaling rules for the time-dependent distribution func- 
tions that characterize the sizes of the infected and the recovered sub-populations. Using heuristic 
arguments, we also obtain scaling laws for the size and duration of the epidemic outbreaks as a 
function of the total population. We perform numerical simulations to verify the scaling predictions 
and discuss the consequences of these scaling laws for near-threshold epidemic outbreaks. 
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I. INTRODUCTION 

The study of epidemics has a long and fascinating his- 
tory that dates back to Daniel Bernoulli who modeled 
the spread of smallpox [3t3|- Theoretical models have 
been quite successful in describing the spread of infec- 
tious diseases. Closely related models have been applied 
to a truly remarkable set of contagious processes includ- 
ing HIV d, [B| , computer viruses @ , spread of technolog- 
ical innovations 0, IH, outbreaks of social and political 
unrest [tj [l(| , and rumor propagation . 

The convenient deterministic framework is most com- 
monly used to model the spread of an epidemic. In this 
formulation, coupled ordinary differential equations de- 
scribe the evolution of macroscopic properties such as 
the average number of infected, the average number of 
recovered, and so on [12]. However, the deterministic ap- 
proach cannot capture fluctuations which are inevitable 
due to the finite size of the population. In many scenarios 
such as the spread of an infection in an isolated boarding 
school with say, a hundred children, statistical fluctua- 
tions are certainly significant. The stochastic framework 
which involves evolution equations for the entire prob- 
ability distribution is instead required to describe finite 
populations [T3-1&]. Theoretical analysis of the stochas- 
tic framework is challenging, even for the most basic in- 
fection processes (l9l-[24| , and many questions concerning 
finite-size effects and extremal properties of the probabil- 
ity distribution functions remain unanswered. 

In this paper, we investigate the stochastic version of 
the classic Susceptible-Infected-Recovered (SIR) infec- 
tion process. In this model, the population consists of 
susceptible, infected, and recovered individuals, and the 
infection spreads through contact between infected and 
susceptible members of the community [12, HH, HE, Hil . 
Each infected individual spreads the infection at a certain 
rate, denoted by a, while infected individuals recover at 
a rate set to one. The epidemic threshold is a = 1. Below 
the threshold, the infection cannot maintain itself; above 



the threshold, the infection can take off. 

We focus on epidemic outbreaks at or near the thresh- 
old [26M32j | - Such "threshold epidemics" are especially 
interesting. From a theoretical viewpoint, the infection 
and the recovery processes are in some sense in perfect 
balance precisely at the threshold. While most infec- 
tions are small, large outbreaks may occasionally happen, 
and hence, threshold epidemics exhibit large fluctuations. 
From a practical viewpoint, human efforts at disease pre- 
vention reduce the infection rate thereby driving infec- 
tious diseases from the pandemic to the endemic phase 
[25|. Meanwhile, natural evolution may increase the in- 
fection rate of diseases hovering just below the thresh- 
old, thereby enhancing the likelihood of an outbreak [33| . 
These epidemics are subtle: they may be difficult to de- 
tect as well as difficult to control. 

We first consider infinite populations and focus on the 
behavior precisely at the epidemic threshold. We ana- 
lyze the rate equation for the probability P^ r (t) that the 
number of infected equals i and the number of recovered 
equals r at time t. The typical number of infected grows 
linearly with time while the typical number of recovered 
grows quadratically with time. We show that the joint 
distribution function obeys the scaling form 



P i>r (t) c±t- 4 $(it-\rt- 2 ) 



(1) 



We obtain the Laplace transform of the scaling function 
$(£, 77) analytically, and present a comprehensive asymp- 
totic analysis. We also discuss scaling properties of the 
respective single-population distribution functions for the 
number of infected or the number of removed. 

Next, we consider finite populations. We combine the 
infinite population results with heuristic arguments to 
derive scaling laws for the size and duration of outbreaks 
in a finite population of size N. At or near the epidemic 
threshold, the effective infection rate is reduced because 
the population is finite. We find that the maximal size 
of the outbreak, M, and the maximal duration of the 
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outbreak, T, obey the nontrivial scaling laws 

M-iV 2/3 , and T~N 1/3 . (2) 

We also formulate the range of infection rates near the 
epidemic threshold for which these scaling laws apply. 

The rest of this paper is organized as follows. We be- 
gin with the definition of the infection process, given in 
Section II. In section III, we treat the infinite popula- 
tion limit by using the rate equation approach. First, we 
derive the outbreak size distribution. We then compute 
two single-variable distributions — the probability dis- 
tribution for the number of infected individuals and the 
probability distribution of the number of recovered indi- 
viduals. Next, we discuss scaling and extremal properties 
of the joint distribution function. In section IV, we con- 
sider finite populations and obtain finite-size scaling laws 
for the size and duration of the outbreaks. We discuss 
the results in section V, and the appendix details some 
necessary inverse Laplace transforms. 



II. THE INFECTION PROCESS 

The spread of infectious diseases has been widely stud- 
ied on regular and irregular spatial structures such as 
lattices [34-42] and complex networks [43M45j ] . Never- 
theless, many characteristics of stochastic epidemics, e.g. 
finite-size properties, remain open questions even for per- 
fectly mixed populations. Throughout this paper we as- 
sume that the populations are well-mixed so that any 
infected individual can spread the disease to any suscep- 
tible individuals. 

In the ubiquitous Susceptible-Infected-Recovered 
(SIR) infection process, the population includes three 
types of of individuals: 

S Uninfected individuals who are susceptible to the 
infection. 

I Infected individuals who are actively spreading the 
disease. 

R Individuals who are neither infected nor susceptible 
including those who have been infected and subse- 
quently recovered, or became immune, or removed. 

An individual may proceed from type S to I to R. In 
this simplified model, all infected individuals may spread 
the disease. 

We investigate the continuous-time version of SIR in- 
fection process [l5j . At any given moment, the popula- 
tion consists of s susceptible, i infected, and r recovered 
individuals; the size of the total population, N = s+i+r, 
remains fixed. The sub-populations change due to two 
competing processes — infection and recovery. An in- 
fected individual may infect a susceptible one with rate 
a/N , where a is the infection rate: 



The overall infection rate is proportional to the size of 
the susceptible population times the size of the infected 
population. Infected individuals recover at a constant 
rate, 



(s,i,r) (s,i- l,r+l). 



(4) 



Without loss of generality, the recovery rate is set to one. 

We consider the natural initial condition where a single 
infected individual is embedded in a "sea" of susceptible 
individuals, that is, (s, i, r) = (N — 1, 1, 0) at time t = 0. 
The infection process ends when no infected individuals 
remain, (s,i,r) — (N — n,0, n). The final number in- 
fected, n, measures the size of the epidemic outbreak. 



III. INFINITE POPULATIONS 

Since the infection process is random, the size of the 
epidemic outbreak is a stochastic quantity. First, we 
investigate the distribution of this quantity as its basic 
characteristics show how the infection process can be in 
one of two phases: an endemic phase where a microscopic 
number of individuals are infected and a pandemic phase 
where a macroscopic number of individuals are infected. 

In this section we analyze the infinite population limit, 
N —> oo. Let G n be the probability that the size of the 
epidemic outbreak equals n. For the initial condition 
(s, i,r) — (N — 1, 1,0), the total infection rate and the 
total recovery rate are a and 1 as follows from ([2])-©. 
In the very first step, the infected individual either recov- 
ers, or a second individual become infected. The former 
happens with probability G% = 1/(1 + a), while the lat- 
ter occurs with the complementary probability a/(l + a). 
Hence, the distribution G n obeys the recursion equation 



G n — 



1 + a 



1 



i-\-j—n 



l 



(•5) 



valid for all n > 1. The convolution term corresponds to 
the case where an additional individual become infected 
— in this situation there are two infection processes, and 
in the infinite population limit these two processes are 
completely independent (4|| 1471 ], Starting with G\ = 
1/(1 + a), the recursion equation gives G 2 = a/(l + a) 3 , 
G 3 = 2a/(l + a) 5 , etc. 

We now introduce the generating function 



n>l 



(6) 



Using the generating function we convert the infinite set 
of recursion equations ([5]) into the quadratic equation, 
(1 + a)Q = aQ 2 + z. Out of the two solutions, the proper 
one satisfies the requirement Q{z = 0) = 0, 



(s,i, 



asi/N 



(s - + l,r). 



(3) 



Qiz) 



1 + a - + a) 2 - 4az 



2a 



(7) 
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The behavior of G(z) in the limit z — > 1 reveals the 
basic characteristics of the infection process. The proba- 
bility Gfinite that only a finite number of individuals are 
infected is given by 



Gfinite — 



1 



a < 1; 
a > 1; 



(8) 



The derivation of (|8]) follows from ([7|) and 
Gfinite = J2n>i G « = ^C 1 )- The threshold value 
a c = 1 separates two regimes of behavior — an endemic 
regime and a pandemic regime. Below the threshold, 
the number of infected individuals is always finite. 
Above the threshold, there is a finite probability that a 
finite fraction of individuals becomes infected, and as a 
consequence, the average size of the epidemic outbreak 
is macroscopic, that is, proportional to the size of the 
total population. 

Also, the average number of infected individu- 
als readily follows from the generating function, viz. 
( rl ) = J2n nG n = £'(!)• Equation (J7J gives 



1 



1 



a < 1. 



(9) 



As expected, the size of the outbreak diverges in the 
vicinity of the threshold. 

To find G n , we expand the generating function G(z) in 
powers of z. The distribution of outbreak size is a prod- 
uct of an exponential and an algebraic factor, expressed 
as a ratio of Gamma functions, 



1 



r n 



(1 



r(n + i) 



1 - 



1 — Q 

1 + a 



By using the asymptotic behavior T(x + a)/T(x) ~ x a 
as x — > oo, we conclude that at the epidemic threshold, 
a = 1 , the size distribution has a power-law tail 



G„^(4^)- 1 /2 n -3/2 ; 



(10) 



as n — > oo. For threshold epidemics, there is balance 
between infection and recovery, as indicated by Q-Q. 
While the majority of outbreaks are small, the algebraic 
behavior (TIT)]) suggests that there is a considerable likeli- 
hood for large outbreaks to occur. 

In the remainder of this section, we focus on the be- 
havior at the epidemic threshold, a = 1. We begin with 
the probability Pi(t) that there are i infected individuals 
at time t. Irrespective of the infection rate, this distri- 
bution function satisfies a closed evolution equation, and 
for the critical a = 1 case, Piit) satisfies 



dt 



(i - l)P;_i + (* + l)P;+i - 2iPi 



(11) 



The initial condition is Pj(0) = Sn. Equations (fTTj) are 
closed because in the infinite population limit, N — > oo, 
the population of infected individuals does not depend on 



the populations of susceptible and recovered individuals. 
We comment that the master equation (fTTj) is a discrete 
diffusion equation with a diffusion coefficient equal to the 
size of the infected population. 

To solve we use the exponential ansatz: 

Pi(t) = *(t) [^{t)} 1 - 1 for i > 0, with the initial condi- 
tions *(0) = 1 and ip(0) = to assure the validity 
of -Pi(O) = Si t ±. This ansatz reduces the infinite set 
of equations (jlip to two ordinary differential equations, 
d^jdt = (1 - ip) 2 and d^/dt = 2(ip - l)ip. Solving these 
coupled equations we obtain 



1 



(1 + *) 2 



t 



l + t 



(12) 



for i > 0. Thus, the size distribution is purely exponen- 
tial. 

The quantity Po(t) gives the probability that the 
infection process has subsided by time t. From 
dPo/dt = Pi and the initial condition Po(0) = 0, we ob- 
tain Po(i) =t/(l + t). Equivalently, one can deduce this 
result from the normalization requirement ^2 i>a Pi = 1. 
Note also that the survival probability P(t), defined as 
the probability that the infection remains active at time 
t, is simply P(t) = 1 — Po(t), and hence, it is given by 



P(t) 



1 



l + t 



(13) 



In the long time limit, the survival probability decays 
algebraically, P ~ t . 

Interestingly, the first moment of the distribution Pi is 
conserved, y\. iPj = 1, and in this sense, the competing 
processes of infection and recovery balance when a = 1 . 
Hence, if we restrict our attention to active outbreaks, 
the average number of infected individuals grows linearly 
with time, (i) = 1 + t, where (i) = J^. iPi/P(t). 

In the long-time limit, the distribution Piit) has the 
scaling form 



Pi(t) ~ t- 2 $(it~ 



(14) 



This scaling behavior immediately follows from (TT21 in 
the limit i — > oo and t — > oo with the scaling variable 
£ = i/t kept fixed. 

The master equation (fTT]> is closed since the infected 
population is decoupled from the other two populations. 
In contrast, the recovered population is coupled to the 
infected population. Therefore, we must analyze the joint 
distribution Pj jT .(t), that is, the probability there are i 
infected and r recovered at time t. Of course, the joint 
distribution function gives a complete description of the 
state of the system and for example, Pi — X)r>o P*,r- The 
joint distribution obeys 



dPj,, 
dt 



(i - l)Pi-i, r + (i + l)P+i,,-i - 2iPi, r (15) 



and Pi : r(0) — 5{ t i d r _Q. In this rate equation, the first gain 
term accounts for infection, while the second gain term 
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represents recovery. For small i and r one can compute 
the joint distribution recursively: Pi o = e~ 2t , Pn 1 = 



(1 - e~ 2t )/2, P M = e- 2t [2t - (1 - e~ 2t )]/2, etc. [46 
Generally, we use the generating function 



(16) 



i>0 r>0 



By multiplying (|15[) by u*w r and summing over all i and 
r, we find that the generating function obeys 



dV 
dt 



= (u 2 -2u + v) 



dV 

du 



(17) 



The initial condition is Vq = V(u,v) = u where 
We can transform equation (|17p into the wave equation 



dV 
dt 



dV 
dw ' 



(18) 



To establish this equation, we introduce the variable 
du 1 



) 2 - 2u' + v 



1 



2yT 



In 



i + vr 



u 1-vT 



1-vT 



« l + Vl 



(19) 



The solution to equation (ITTl) is simply Vo{w + t). Since 
Po = u, we need to obtain u in terms of the variable w 
and then, replace w with u; + t to obtain the generating 
function explicitly. From (j!9[) . the expression 



= 1 - VI - u 



2^ 



1 + yT^y 
l - vT"^' 



(20) 



gives Po = u i n terms of v and u>. As a function of the 
variables w and v, the generating function becomes 



V(w,v) = 1 - VT 



2Vl _r ~ 



1 + e 2(«,+i)VT=T 
1 - VT - v 



-, -l 



1 



This result is derived from (|20l) by replacing w with lo+f. 
In terms of the original variables, the generating function 
is 



V(u,v) = 1 - vT 



2VT 



l-VT 



D-(E-l)u 



(21) 



We obtained this result by expressing exp fewyi — v\ in 
terms of u. In writing (|2Tj) . we also used the shorthand 
notations 



E = e z 



D={E+l)VT zr v + E-l. (22) 



For example, we can verify that the generating function 
yields the size distribution of outbreaks. In the long- 
time limit, the last term on the right-hand-side of (f2"Tj) 



vanishes since the quantities D and E are both divergent. 
With the shorthand notation P^ = limt-^x, P, we have 
Voo(u,v) = 1 — VI ~ v - By substituting a = 1 into 0, 
we confirm that Pi, r (t) — > 5j ; o G r when i — > oo. 

First, we analyze the probability distribution Tl r (t). 
By definition, LI r (t) is the probability to have r recovered 
individuals at time t. Of course, II r (i) = ^2i> Pi,r(t). 
The distribution of recovered differs from the distribution 
of infected in that in the long-time limit, it approaches a 
nonzero value, II r — > G r , whereas Pi — > for all i > 0. 
To study the long-time asymptotic behavior, we write 
II r (t) = G r + H r {t). We have H r (t) -> as t -> oo. The 
corresponding generating function -ff(u) = X)r>o ^r" r is 
given by 



JJ(«) 



2VT^=U 



1 



(23) 



We obtain this expression from the joint generating func- 
tion (jnj by setting u = 1, H (v) = V(l,v) - Poo(l,«) 
where P<x>(l, u) = 1 — \A - V. 

According to the scaling behavior (fT4"|) , the typical size 
of the infected population grows linearly with time, i ~ t. 
Heuristically, dr/dt ~ z since the recovery rate is con- 
stant. Consequently, the typical size of the recovered 
population is quadratic, r ~ i 2 , and hence, we expect 
the scaling behavior 



H r {t) ~ t~ 3 cp(r i" 2 ). 



(24) 



The Laplace transform of the scaling function ip(rj) fol- 
lows immediately from behavior of the generating func- 
tion (f23|) in the limit v — >• 1. We take the limits v —¥ \ 
and i — > oo while keeping the variable b — t 2 (l — v) fixed. 
In this limit, we have v r — > e~ bv and the sum over r turns 
into the integral 



r>0 



dr\ e 



- bi] 



By using this scaling transformation along with Eq 
we find the Laplace transform of the scaling function 



drje br ><p(ri) = ^ 
o e 2 ^ + 1 



(25) 



The inversion of this Laplace transform through inte- 
gration in the complex plane is detailed in the Appendix 
where we show that the scaling function ip(r)) is given by 



^)^2^£(fc + I) 2 e - 2 (^) 2 "--i 



k=0 



(26) 



We now substitute this expression into the scaling form 
(f24f and observe that the algebraic factor (47r) _1 / 2 77 _3/ ' 2 
and the final distribution G r cancel each other. Thus, 
the distribution H r (t) has the same scaling form as the 
distribution H r (t), 



IL,(i) ~t~ z cj>(rt- 2 ). 



(27) 
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The scaling function cj>(j)) is given by the sum 

oo 

<M77) = 2^^(fc+i) 2 e -* 2 ( fe +5)\ (28) 



k=0 



The leading asymptotic behaviors of the scaling function 
4>(rf) are 

f(4 7 r)-V2 7? -3/ 2 

W) - W_2 /9V-tV4 ( 2 °) 



(^/2)e 



?7 — > oo. 



The algebraic behavior in the small-77 limit is consistent 
with the power-law tail (jlOp . To obtain this behavior, we 
simply convert the sum in (|28[) into an integral. We note 
that with the scaling form (l27l) and the small-?? diver- 
gence, the quantity ^ r II r is indeed finite. Also, the first 
term in the series yields the exponential behavior in the 
large- 77 limit. Hence, both the distribution of recovered 
and the distribution of infected have exponential tails. 

We now analyze the joint distribution Pi >r {t). We re- 
strict our attention to active epidemics that correspond 
to i > 0. Thus, we focus on the following component of 
the joint generating function 



U V . 



(30) 



This component follows directly from the generating 
function, T(u,v) — V(u,v) — V(0,v), and by using the 
explicit solution (j2"Tj) . we obtain 



J-"(u, v) = 



AEu(l - v) 
D[D — (E-l)u]' 



(31) 



The scaling behaviors (|14l) and (|2"T|) for the single- 
population distributions imply that the joint distribu- 
tion function adheres to the scaling form ([1]) stated in 
the introduction. Since the survival probability P{t) = 
J2i>i J2 r Pi;r(t) decays with time according to ([TB")) . the 
scaling function is normalized J J dr/dt; $(£, 77) = 1. 

The above analysis suggests use of the joint Laplace 
transform 



b)= / dtdrie-^Q&r,), (32) 
Jo Jo 



that is merely the continuous counterpart of the gener- 
ating function. To obtain f(a,b) from F(u,v) given in 
(1311) . we take the limits t — ¥ 00, u — > 1, and v —¥ 1 while 
keeping the variables 

a = (1 — u) t and b = (1 — v) t 2 

fixed. By taking these limits, we observe that the right- 
hand side of (|30|) becomes t~ 1 f{a, b), and find the joint 
Laplace transform in explicit form 



f(a,b) 



Vb 



1 



sinh vb 



Vb coth Vb 



(33) 



We stress that the quantity /(a, 6) describes only active 
infection processes. One can verify that /(0, 0) = 1. 
Moreover, in the b — > limit, we have /(a, 0) = 1/(1 + a), 
for which the inverse Laplace transform is immediate, 
J drj $(£,77) = e - ^. Indeed, we recover the scaling func- 
tion $(£) in Eq. (fT4|). 

Since the inverse Laplace transform of /(a, b) with re- 
spect to the variable b is immediate, we have 



drje- bri ^, V ) 



Vb 



sinh Vb , 



-£\/6coth \/b 



(34) 



The Appendix outlines how to invert this Laplace trans- 
form to obtain the leading asymptotic behaviors for large- 
77 and small-77, 



4(l+g/2) 3 

2 ,1/4 



exp 



(l+g/2) 2 
V 



0. 



exp [— 7r 2 ?7 + n\/8^rj — £] rj — > 00. 

These limiting behaviors apply for a fixed value of £. 

The function /(a, 6) captures all moments of the scal- 
ing function $(£, rf). For example, by expanding equation 
([33| as a Taylor series in the conjugate variables a and 
6, /(a, 6) = 1 — a — |fo + ab + ■ ■ ■ , and comparing with 
the definition (I32[) . we obtain the lowest-order moments 
(£) = 1, (77) = I, and (£77) = 1. In particular, 



(36) 



so there is positive correlation between the size of the 
two populations. Intuitively, we expect that long-lasting 
epidemic outbreaks involve large numbers of infected and 
recovered, while the opposite is true for short-lived out- 
breaks. 

For completeness, we mention that the above anal- 
ysis can be repeated for inactive outbreaks. Starting 
with "P(0, v) given in ([2T]) and following the steps lead- 
ing to (|2"7)) . we find that the distribution Po,r(i) that 
an epidemic outbreak has ended by time t and that 
the size of the outbreak equals r, has the scaling form 
Po,r(i) — t~ 3 4>(r t~ 2 ). The scaling function ^(77) resem- 
bles 4>(r]) given in 



</>(77) = 2tt 2 fc2 e 



2 1 2 
-7T K TJ 



(37) 



k=l 



The leading asymptotic behavior in the small-77 limit 
is identical to that in (|2"9|) . and again, there is 
exponential decay, albeit with twice the coefficient, 
(j>(r)) ~ 2-7T 2 exp(— 7r 2 ?7) in the large-77 limit. As expected, 
inactive outbreaks dominate the size distribution of the 
recovered population at large times. 



IV. FINITE POPULATIONS 

The zeroth and first-order moments of the distribution 
G r given in equations © and (JSJ) show that the size of 
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FIG. 1: The average outbreak size versus the population 
size for the SIR infection process at the epidemic threshold 
(a = 1). Shown are Monte Carlo simulation results repre- 
senting an average over 10 9 independent realizations of the 
infection process (circles). A line of slope 1/3 is also shown 
as a reference. 
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FIG. 2: The average outbreak duration at the epidemic 
threshold versus the population size. Simulation results, ob- 
tained from an average over 10 9 realizations (circles) are com- 
pared with the reference line (1/3) In N + C. 



the epidemic outbreak is microscopic below the epidemic 
threshold, a < 1, but macroscopic above the threshold, 
a > 1. These results hold for infinite populations, yet 
real- world epidemic outbreaks involve finite populations, 
and there is practical need for understanding how basic 
characteristics such as the average size and the average 
duration of the epidemic depend on the size of the pop- 
ulation, N. 

Finite-size effects are most pronounced in the vicinity 
of the epidemic threshold. Let us consider large but finite 
populations, N ^> 1, and let us consider a scenario in 
which the susceptible population has been depleted by 
n, that is, s — N — n. From the very definition of the 
infection process ([3]) , we determine that the infection rate 
is reduced, a — > a*(iV), with 

a*(JV)=a(l--|), (38) 

because the total population is finite. Clearly, as the 
large "reservoir" of susceptible individuals shrinks, the 
infection process weakens. 

Let us now consider the average size of the outbreak for 
a threshold epidemic (a = 1). We assume that there is 
a maximal outbreak size M, and that the outbreak can 
not exceed this size due to depletion in the number of 
susceptible individuals. On the one hand, equations (|38[) 
and © suggest that (n) ~ N/M. On the other hand, 
the algebraic distribution (fTU)) gives a second estimate 
(n) ~ X) n <M n_1/2 ~ M 1 / 2 for the average size of the 
outbreak. Equating these two estimates, N/M ~ M 1 ' 2 , 
we conclude that M ~ TV 2 / 3 as stated in (2). While a 
naive interpretation of the effective infection rate (|38|) 
would suggest that the cutoff is proportional to the to- 
tal population, we find that threshold epidemics have a 
substantially smaller upper bound. 



Using (n) ~ N/M and M ~ N 2 ' 3 we see that the 
average size of a threshold epidemic grows as 

(n) ~ iV 1 / 3 . (39) 

This interesting scaling law was established by Martin- 
Lof fl9ll and it has been confirmed in several recent stud- 
ies |20h24|] . Figure [1] shows a numerical verification of 
this behavior. For finite populations, this scaling law 
holds in the neighborhood (often termed critical or scal- 
ing windows) of the epidemic threshold. To estimate the 
size of this neighborhood, we simply compare (|3T)f with 
©. The size of the outbreak grows sub-linearly in the 
population size as long as 

|l-a| ~A^ 1/3 . (40) 

Of course, this neighborhood shrinks as the size of the 
population grows. Yet, for moderate populations, the 
size of this neighborhood is considerable. 

The scaling law (f3T?| indicates that for finite popula- 
tions, there are actually three regimes of behavior. Well 
below the epidemic threshold, a finite number of individ- 
uals becomes infected. Well above the threshold, a finite 
fraction of the population becomes infected. In the vicin- 
ity of the threshold, the size of the outbreak grows as the 
1/3 power of the population size, 

(0(1) C^oo, 
(n>~ InV*Y(0 C = 0(1), (41) 

[N C^-oo, 

with the scaling variable ( = iV 1 ^ 3 (l — a). A finite value 
of C indicates that the infection rate is in the neighbor- 
hood of the threshold. (The scaling function Y(() can be 
extracted from Refs. [lj, HjlH.) 
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FIG. 3: The survival probability at the epidemic threshold. 
Shown is the normalized survival probability P(t,N) / P(t) 
versus the normalized duration time t/N 1 / 3 . The data corre- 
sponds to an average over 10 8 realizations. 



The duration of a threshold epidemic also adheres to a 
finite-size scaling law. Since the size of the recovered pop- 
ulation grows quadratically with time, r ~ t 2 , as implied 
for example by the scaling behavior (|14[) . we conclude 
that there is a maximal time scale T, and that the du- 
ration of a near-threshold outbreak may not exceed this 
scale. We can estimate the time scale T ~ N 1 / 3 stated in 
@, from the heuristic argument T 2 ~ M. The average 
duration of the infection process follows from the survival 
probability, (t) = — f Q dttdP/dt, and using P ~ i -1 , we 
find that the average duration grows logarithmically with 
time (Fig. 12), 



(t) ~ - In AT. 

v 7 3 



(42) 



Logarithmic growth, albeit with a unit prefactor, was 
predicted by Ridler-Rowe (iij . 

The scales M ~ TV 2 / 3 and T - N 1 / 3 fully character- 
ize the size and duration of the infection process when 
a = 1. For example, the survival probability P(t,N) to 
have at least one infected at time t obeys P(t, N)/P(t) — > 
S (t/N 1 / 3 ) (Fig. [3J . However, the simulations show that 
the convergence to this scaling form is not uniform — it 
is slow for short durations but fast at large durations. 

We performed numerical simulations to check the scal- 
ing predictions for the case a = 1. Our numerical method 
is merely a Monte Carlo procedure to solve the mas- 
ter equation (fT5|) . We conveniently keep track of two 
populations, s and i, from which the overall infection 
rate and recovery rate are respectively Ri = si/N and 
R r = i. The probabilities Pj = Ri/(Ri + R r ) and 
P r = R r /(Ri + R r ) that the following step involves infec- 
tion or recovery are calculated and then, the populations 



(«!*) -> 



(*-l,i + l) 
(M-l) 



with probability Pi, 
with probability P r . 



(43) 



The infection process ends when i = for the first 
time. Time is updated by the inverse of the total rate 
t — > t + l/(Ri + Rr)- With this straightforward numeri- 
cal procedure, we can simulate as many as 10 9 indepen- 
dent runs in populations as large as N = 10 9 . As shown 
in figures Q] and [2] there is excellent agreement between 
the numerical results and the theoretical predictions for 
the average size (I3U1) and the average duration (|4*2"j) . 



V. DISCUSSION 

We investigated various statistical properties of the 
SIR infection process at a threshold infection rate. We 
analyzed the rate equations for the two-population dis- 
tribution function that characterizes the probability that 
the system has a specified number of infected and recov- 
ered individuals. Our analysis yields scaling behavior in 
the asymptotic long-time limit and gives extremal prop- 
erties of the joint distribution function as well as the as- 
sociated single-population distributions. We used these 
infinite population results to justify scaling laws for finite 
populations. 

Outbreaks in the vicinity of the epidemic threshold 
have a distinct size, characterized by non-trivial power- 
law dependence on the total population size. This scaling 
behavior applies in the vicinity of the epidemic thresh- 
old. The size of this neighborhood, N" 1 / 3 , is larger than 
the canonical value iV -1 / 2 expected from the traditional 
large-size expansions or from the deterministic descrip- 
tion Therefore, statistical fluctuations and finite 
population effects are pronounced and may be subtle 
near the epidemic threshold. We note the an identical 
scaling arises near the percolation point for Erdos-Renyi 
random graphs [SJ [52[ • Additional connections between 
SIR infection processes and random graphs were reported 
recently [2a |. 

The scaling laws for the time dependence and the size 
dependence are useful. For example, the scaling laws 
for the critical kinetics can be used to find the epi- 
demic threshold in numerical simulations of infection pro- 
cesses on complex networks for which the threshold is not 
known analytically. Furthermore, the number of coupled 
ordinary equations needed to compute the joint distri- 
bution Pi tr numerically is in principle quadratic with N . 
However, the scaling laws i ~ N 1 / 3 and r ~ N 2 / 3 imply 
that the relevant number of equations is much smaller, 
being proportional to N. 

For large but finite populations, we understand the ba- 
sic scaling laws, but much less is known about finite-size 
scaling functions. The only exception is the known scal- 
ing function that gives the size distribution of the out- 
breaks at the epidemic threshold [l9l . l2~i| . The analytic 



8 



determination of the scaling function characterizing the 
duration of outbreaks near the epidemic threshold is a 
challenging problem because there is no closed equation 
for the total duration of the outbreak [l8[ . 

Finally, we mention that in this study we assumed 
that all individuals can interact. In most applications, 
the spatial [3^ - l42j or network [Zij - flsj structure of the 
infected domain play an important role. Finding the 
corresponding scaling functions for the time-dependent 
behavior at the critical point or for the finite-size scaling 



are also challenging open problems. 
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Appendix A: Inverse Laplace Transforms 

From Eq. (|25|) , the scaling function ip(r]) equals the 
inverse Laplace transform 



7+ioo 



dbe h ^- 



2Vb 



7—200 



aVb 



+ 1 



(Al) 



where the integration contour is a line parallel to the 
imaginary axis in the complex b plane and 7 > so that 
all singularities are to the left of the integration contour. 
The integrand has simple poles at bk = —tt 2 (k + i) , 
k = 0,1,2,... and a branch point at the origin, and we 
select a branch cut along the negative real axis so that it 
does not cross the path of integration. Furthermore, we 
pick a closed Bromwich contour formed by the contour 
followed by a quarter of a circle of infinitely large ra- 
dius, followed by the top of the branch cut with infinites- 
imal half-circles around each pole and then encircling the 
origin and proceeding similarly along the bottom of the 
branch cut, followed by a quarter of a circle. By the 
Cauchy theorem, the integral along this closed contour 
vanishes. The integrals over the circles vanish and the 
integrals over the branch cut can be computed to yield 



^) = 27T 2 V(fc+i) 



fc=0 



,2 -7 

e 



- I dce- CT >- 



where the sum is proportional to the residues of poles at 



b k 



Since the integral on the right-hand 



side equals (47r) _1 / 2 ?7 _3 / 2 , we arrive at (|26l) . 

The inverse Laplace transform of (|34j) with respect to 
b gives the joint scaling function 



1 ry+ico 

F &ri) = — dbe b ^ 

£Tt1 ./ 7 _ioo 



Vb 



sinh y/b 



' Vb coth Vb 



To obtain the leading asymptotic behavior in the small-?? 
limit, given in (1351) . we simply evaluate the large- b be- 
havior. In the opposite limit, 77 — > 00, the asymptotic be- 
havior follows from the singularity of the Laplace trans- 
form in the complex b plane that is closest to the origin. 
The Laplace transform has singularities at bk = —ir 2 k 2 , 
k = 1,2,..., so we focus on the singularity at b\ = — 7r 2 . 
For a fixed £, the scaling function appears to decay ex- 
ponentially, 77) ~ e _7r n as 77 — > 00. To establish the 
complete asymptotic behavior we recall that a subdom- 
inant power-law prefactor, F(£,r)) ~ rj m e _7r n , implies 
that the Laplace transform has the algebraic singularity 
(•7T 2 + 6)~ m ~ 1 as b —> —tt 2 . In fact, the Laplace transform 
has the essential singularity (we set b — — ir 2 (l — e) 2 , so 
that the limit b — > —tt 2 is equivalent to the e — > limit) 



e" 2 exp [(e- 1 - 



(A2) 



An essential singularity of the type e^ e corresponds to 
a subdominant prefactor that is a stretched exponential, 
e nV8f;?7_ A n additional power-law factor times a numer- 
ical factor allow to match the singularity (| A2|) . The re- 
sulting large- 77 asymptotic behavior is given in (|35[) . 



